1 Carga de paquetes

for (package in c("tidyverse","fpp3", "GGally", "normtest")) {
    if (!require(package, character.only=T, quietly=T)) {
        install.packages(package)
        library(package, character.only=T)
    }
}
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## ── Attaching packages ──────────────────────────────────────────── fpp3 0.4.0 ──
## ✓ lubridate   1.8.0     ✓ feasts      0.2.2
## ✓ tsibble     1.1.1     ✓ fable       0.3.1
## ✓ tsibbledata 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date()    masks base::date()
## x dplyr::filter()      masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval()  masks lubridate::interval()
## x dplyr::lag()         masks stats::lag()
## x tsibble::setdiff()   masks base::setdiff()
## x tsibble::union()     masks base::union()
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

2 Funciones auxiliares

# Contraste para los coeficientes
my_t_test <- function (object, ...) 
{
  par <- rbind(t_stat=tidy(object)$statistic, p_value=tidy(object)$p.value)
  colnames(par) <- tidy(object)$term
  if (NCOL(par) > 0) {
    cat("\nt-test:\n")
    coef <- round(par, digits = 4)
    print.default(coef, print.gap = 2)
  }
}

# Gráfico de correlogramas de residuos
my_tsresiduals <- function (data, ...) {
  if (!fabletools::is_mable(data)) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  data <- stats::residuals(data)
  if (n_keys(data) > 1) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial",  
    ...)
}

# Validación cruzada anidada

nested_cv <- function(df, h, last_train, string_formula){
  
  nested_errors <- vector()  
  
for (i in seq(last_train, last(df$date), h)){
  train <- df %>% filter(date<=i)
  test <- df %>% filter(date>i)
  
fitted_model <- train %>%
  model(arima = ARIMA(as.formula(string_formula)))

h_forecast = min(dim(test)[1], h)

fc <- fitted_model %>%
  forecast(h=h_forecast)

test_err <- fc %>%
  accuracy(test) %>%
  select(MAPE)
nested_errors <- c(nested_errors, test_err$MAPE)
NewList <- list("errors"=nested_errors, "mean"=mean(nested_errors))
}
return(NewList)}

# Test de autocorrelacines de los residuos
autocorrelation_test_plot <- function(aug, dof = 4, m = 7,  h = 5, alpha = 0.05){
vec <- c()
num_lags = seq(1, h*m)

for (i in num_lags){
 vec <- c(vec,aug %>% features(.resid, ljung_box, lag=i, dof=dof) %>% .$lb_pvalue)
}

autocorr_pvalues_resid <- tibble(
  lag = num_lags, 
  p_value = vec,
  incorelated = p_value >= alpha
)

plot <- autocorr_pvalues_resid %>% 
  drop_na() %>% 
   ggplot(aes(lag, p_value, color = incorelated)) + 
  geom_point() + 
  geom_hline(aes(yintercept = alpha), linetype="dashed", color = "indianred2")


newList <- list("values" = autocorr_pvalues_resid, "plot" = plot)
return(newList)
}

3 División train y test

Antes de empezar a hacer el esco2dio hay que distinguir la parte de los datos que se utilizará para construir la fórmula (para modelizar) de aquella que se utilizará para validar los resultados y cuyos datos, no deben ser utilizados.

Convertimos la serie en un objeto tsibble, para trabajar de manera cómoda. Analizamos visualmente la serie y las marcas de tiempo para saber cómo dividir los datos.

co2 <- read_csv('co2_mm_mlo.csv')
## Rows: 746 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): year, month, value
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
co2$date <- with(co2, sprintf("%d-%02d", year, month)) # juntar mes y año
co2$date = tsibble::yearmonth(co2$date)
co2 = as_tsibble(co2, index = date)


co2 %>%
  autoplot(value) +
    labs(title = "CO2") +
    xlab("Fecha") + ylab("   ") +
  geom_point(aes(y = value), size = 1.2, shape = 16, color = "black")

Dado que tenemos datos mensuales y visualmente se aprecia una estacionalidad anual, reservamos el último año para test y el resto para train.

co2_train <- co2 %>% filter_index(. ~ "2005-12-01")
co2_test <- co2 %>% filter_index("2006-01-01" ~ .)

4 Fase de identificación

Se comienza la fase de identificación donde se deciden las transformaciones a realizar y el propio ajuste de la serie. Aunque veremos que cuando se realiza la selección de los parámetros del modelo ARIMA, la iteración con la estimación y el contraste son constantes hasta converger a un modelo válido.

4.1 Estacionariedad

Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada.

co2_train %>%
gg_season(value, period = "year",labels = "left")

Estacionalidad anual: la forma de las series anuales son iguales.

Descomposición STL para analizar la existencia de tendencia y estacionalidad.

co2_train %>%
   model(seats = feasts:::STL(value)) %>%
  components() %>%
  autoplot()

Se confirman visualmente la estacionalidad anual y una tendencia creciente y luego decreciente.

A continuación, esco2diamos si la serie es estacionaria en media y varianza y acco2amos en consecuencia.

Estacionariedad en varianza

Analicemos qué transformación Box-Cox puede aplicarse para estabilizar la varianza.

lambda <- co2_train %>%
features(value, features = guerrero) %>% pull(lambda_guerrero)
lambda
## [1] 0.7022992
co2_train %>% autoplot(box_cox(value,lambda)) +
labs(y = "Box-Cox transformed TU")

Dado que la serie visualmente muestra una varianza más o menos constante, que el valor Box Cox es más próximo a 1 que a otros valores de transformaciones usuales (logaritmo, raíz cuadrada) y que no ha habido un cambio muy llamativo con la transformación, no realizaremos ningún cambio a priori. Si en el proceso de estimación no conseguimos un modelo adecuado, volveríamos a este punto para valorar de nuevo la transformación.

Estacionariedad en media

La serie muestra una tendencia creciente y luego decreciente, es decir, la media va variando a lo largo de la serie. Por tanto, es esperable que se necesite realizar al menos una diferencia para estabilizar la media. Formalmente, se comprueba mediante una prueba de raíces unitarias.

co2_train %>%
  features(value, unitroot_kpss)

El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada. Es decir, hay evidencias estadísticas para rechazar la hipótesis nula de que la serie sea estacionaria. Por tanto, habría que realizar una diferencia regular.

co2_train %>%
  features(difference(value, 1), unitroot_kpss)

Una vez realizada la diferencia regular, el p-valor del test indica que no hay evidencias estadísticas para rechazar la hipótesis de que la serie sea estacionaria en media.

Graficamos la serie diferenciada y ya observamos que el nivel de la serie no cambia:

co2_train %>% autoplot(difference(value, 1)) + labs(title = "Serie con una diferencia regular")

Debido a la estacionalidad de los datos, quizá deba realizarse una diferencia estacional. Para evaluarlo podemos recurrir a la función unitroot_nsdiffs(), esta evalúa el estadístico fuerza estacional \(F_S\) y sugiere una diferencia estacional si esta es mayor que 0.64.

co2_train %>%
  mutate(turnover = difference(value,1)) %>%
  features(turnover, unitroot_nsdiffs)

El método sugiere una diferencia estacional. En cualquier caso, para asegurar que la serie es estacionaria, en la fase de identificación del modelo debemos asegurarnos de que las raíces de la parte autoregresiva están fuera del círculo unitario (y, por tanto, 1/raíz dentro del círculo).

4.2 Determinación del modelo

En los pasos previos se ha sugerido una diferencia regular y una estacional. Empezaremos teniendo en cuenta con la diferencia regular. Es decir, nos dirigimos a identificar los órdenes \(p, q, P, Q\) de un SARIMA (p,d,q)x(P,D,Q)\(_{12}\) donde ya sabemos que \(d=1\). Vemos que los residuos contienen mucha información.

fit <- co2_train %>%
  model(arima = ARIMA(value ~  1+pdq(0,0,0) + PDQ(0,0,0, period = 12)))
fit %>% my_tsresiduals(lag_max = 36)

Claramente \(p=1\).

5 Fase de estimación y contraste

Comencemos ajustando la parte estacional mediante un SARIMA (1,0,0)x(0,0,0)\(_{12}\).

fit2 <- co2_train %>%
  model(arima = ARIMA(value ~ 1+pdq(1,0,0) + PDQ(0,0,0, period=12)))
fit2 %>% my_tsresiduals(lag_max =36)

report(fit2)
## Series: value 
## Model: ARIMA(1,0,0) w/ mean 
## 
## Coefficients:
##          ar1  constant
##       0.9991    0.3188
## s.e.  0.0012    0.0232
## 
## sigma^2 estimated as 1.479:  log likelihood=-893.43
## AIC=1792.86   AICc=1792.91   BIC=1805.8
my_t_test(fit2)
## 
## t-test:
##               ar1  constant
## t_stat   853.9515   13.7145
## p_value    0.0000    0.0000

Este valor indica que hay que diferenciar (como ya habíamos comprobado).

Diferenciamos, \(d=1\). SARIMA (0,1,0)x(0,0,0)\(_{12}\).

fit3 <- co2_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(0,0,0, period=12)))
fit3 %>% my_tsresiduals(lag_max =36)

report(fit3)
## Series: value 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         0.1159
## s.e.    0.0515
## 
## sigma^2 estimated as 1.463:  log likelihood=-886.2
## AIC=1776.4   AICc=1776.43   BIC=1785.03
my_t_test(fit3)
## 
## t-test:
##          constant
## t_stat     2.2518
## p_value    0.0247

La ACF muestra una estructura AR estacional \(P=1\). Es decir, SARIMA (0,1,0)x(1,0,0)\(_{12}\).

fit4 <- co2_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(1,0,0, period=12)))
fit4 %>% my_tsresiduals(lag_max =36)

report(fit4)
## Series: value 
## Model: ARIMA(0,1,0)(1,0,0)[12] w/ drift 
## 
## Coefficients:
##         sar1  constant
##       0.9436    0.0068
## s.e.  0.0124    0.0131
## 
## sigma^2 estimated as 0.1652:  log likelihood=-297.76
## AIC=601.52   AICc=601.56   BIC=614.45
my_t_test(fit4)
## 
## t-test:
##             sar1  constant
## t_stat   76.1768    0.5210
## p_value   0.0000    0.6026

Introducimos \(Q=1\) (se salen los múltiplos de 12 en la PACF y la correlación más alta en la ACF es la de 12). Es necesario también hacer una diferencia estacional.

SARIMA (0,1,0)x(0,1,1)\(_{12}\).

fit5 <- co2_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(0,1,1, period=12)))
fit5 %>% my_tsresiduals(lag_max =36)

report(fit5)
## Series: value 
## Model: ARIMA(0,1,0)(0,1,1)[12] w/ poly 
## 
## Coefficients:
##          sma1  constant
##       -0.9200    0.0024
## s.e.   0.0204    0.0015
## 
## sigma^2 estimated as 0.09687:  log likelihood=-139.1
## AIC=284.2   AICc=284.24   BIC=297.06
my_t_test(fit5)
## 
## t-test:
##             sma1  constant
## t_stat   -45.108    1.6117
## p_value    0.000    0.1076

Habría que meter \(p=1\), \(q=1\). SARIMA (0,1,1)x(0,1,1)\(_{12}\).

fit6 <- co2_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(0,1,1, period=12)))
fit6 %>% my_tsresiduals(lag_max =36)

report(fit6)
## Series: value 
## Model: ARIMA(0,1,1)(0,1,1)[12] w/ poly 
## 
## Coefficients:
##          ma1     sma1  constant
##       -0.346  -0.8874    0.0024
## s.e.   0.044   0.0216    0.0012
## 
## sigma^2 estimated as 0.08865:  log likelihood=-112.16
## AIC=232.33   AICc=232.4   BIC=249.49
my_t_test(fit6)
## 
## t-test:
##              ma1      sma1  constant
## t_stat   -7.8643  -41.1472    2.0855
## p_value   0.0000    0.0000    0.0375
gg_arma(fit6)

5.1 Diagnosis de residuos

Para completar la fase de contraste evaluamos los residuos. Veamos primero el histograma.

aug <-fit6 %>% augment()

# Histogram
aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram(bins = 50) +
  ggtitle("Histogram of residuals")

Test media = 0

Contraste t-student para contrastar si la media es 0.

# Student's t-Test for mean=0
t.test(aug$.resid)
## 
##  One Sample t-test
## 
## data:  aug$.resid
## t = 0.0027261, df = 551, p-value = 0.9978
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.02451750  0.02458564
## sample estimates:
##    mean of x 
## 3.407357e-05

Como p-value > 0.05 no existen evidencias significativas para rechazar la hipótesis nula de que la muestra tenga media 0.

Test autocorrelaciones

A continuación comprobamos que los residuos están incorrelados.

# Ljung-Box autocorrelation
aug %>% features(.resid, ljung_box, lag=12, dof=3)

Podemos concluir que los residuos están incorrelados.

Test homocedasticidad

Para contrastar la heterocedasticidad se puede utilizar una regresión media-dispersión. Calculando los grupos de manera anual, dado que esa es la estacionalidad de la serie.

log_log <- aug %>% as_tibble() %>% 
  group_by(year(date)) %>% 
  summarize(mean_resid = log(mean(.resid+1)), std_resid = log(sd(.resid+1))) 
  
  summary(lm(std_resid~mean_resid, log_log))
## 
## Call:
## lm(formula = std_resid ~ mean_resid, data = log_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.25683 -0.14625  0.01988  0.20569  0.39826 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.27734    0.04477 -28.530   <2e-16 ***
## mean_resid  -0.78133    0.71330  -1.095    0.279    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3035 on 44 degrees of freedom
## Multiple R-squared:  0.02655,    Adjusted R-squared:  0.004422 
## F-statistic:   1.2 on 1 and 44 DF,  p-value: 0.2793

El p-valor de log(media) es mayor que 0.05, no hay evidencias para rechazar la hipótesis de que sean homocedásticos.

Test de normalidad

Para contrastar la normalidad realizamos el test Jarque-Bera para evaluar la normalidad.

# Jarque Bera test
jb.norm.test(na.omit(aug$.resid))
## 
##  Jarque-Bera test for normality
## 
## data:  na.omit(aug$.resid)
## JB = 9.0477, p-value = 0.016

Como el p-valor <0.05 se rechaza la hipótesis nula de normalidad de los residuos.

6 Fase de predicción

Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones en el conjunto de test, calculamos sus errores y los comparamos con los residuos.

# Residual accuracy
resids <- fit6 %>% 
  accuracy() %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Training') 

# Forecasting
fc <- fit6 %>%
  forecast(h=12) 

test_err <- fc %>% 
  accuracy(co2) %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Test')

# Show errors together
bind_rows(test_err, resids) %>% select(Evaluation, everything())

Además de evaluar los errores con estas medidas globales, podemos evaluar los errores cometidos año a año visualmente. Así, podemos detectar outliers o efectos de calendario no detectados. Pudiendo así corregirlos con modelos más complejos mediante variables exógenas.

aug %>%  ggplot() +
  geom_line(aes(x = date, y = .fitted), color="red") +
  geom_line(aes(x = date, y = value), color="gray24") +
  ggtitle("SARIMA train fitted values") +
  xlab('Dates') +
  ylab(' ') #+ facet_wrap(vars(year(date)), scales = 'free')

También podemos dibujar las predicciones con los intervalos de confianza (aunque en este caso no sean representativos, dada la no normalidad de los residuos). Evaluamos así también la capacidad de generalización del modelo, comparando las predicciones con el conjunto de test.

h <- dim(co2_test)[1]

demanda_plot <- co2 %>% filter(date>last(co2_train$date)-14 & date< last(co2_train$date) + 120 )

fit6 %>%
  forecast(h=120) %>%
  autoplot(demanda_plot)

Como sabemos que los modelos SARIMA no predicen bien más alla del orden más alto disponible, en este caso 12, los errores en el conjunto de test más allá de 12 datos se dipararían.